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We investigate energy transport in several two-level atom or spin- 1/2 models by a direct coupling 
to heat baths of different temperatures. The analysis is carried out on the basis of a recently derived 
quantum master equation which describes the nonequilibrium properties of internally weakly coupled 
systems appropriately. For the computation of the stationary state of the dynamical equations, we 
employ a Monte Carlo wave-function approach. The analysis directly indicates normal diffusive or 
ballistic transport in finite models and hints toward an extrapolation of the transport behavior of 
infinite models. 
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I. INTRODUCTION 

The transport of energy or heat has been intensively 
studied since Fourier introduced his famous law of heat 
conduction in 1807. Surprisingly, still 200 years later, 
some fundamental problems remain unsolved- . Contrary 
to our everyday experience, the appearance of diffusive 
behavior according to Fourier's famous law is difficult to 
obtain from the direction of any underlying microscopic 
theory. 

A question of central relevance concerns the classifi- 
cation of the transport properties of a system into nor- 
mal diffusive or ballistic behavior. In the classical do- 
main, it seems to be largely accepted that normal en- 
ergy transport, i.e., spatial diffusion instead of ballistic 
transport or localization, requires the chaotic dynamics 
of a nonintegrable system 2 -^. In the quantum regime, 
however, the question whether diffusive behavior follows 
from the underlying theory turns out to be a controver- 
sial issue^iSiZiS. This is mostly due to the nontrivial 
character of the question, how energy is transported on 
the microscopic scale. 

Amongst the many different techniques of investigating 
transport in quantum mechanics, let us consider two ap- 
proaches in more detail here. The first one is the promi- 
nent Green-Kubo formula. Derived on the basis of linear 
response theory it has originally been formulated for elec- 
trical transport ^ 10 ' 11 . Therein, the system is perturbed 
by an external force, first the electric field, applied to 
the system. The resulting current of charge is viewed 
as the response to this external perturbation. Finally, 
the transport coefficient (conductivity) follows from a 
current-current autocorrelation. The same approach is 
also used in the different case of density driven transport, 
e.g., the transport of energy or heat. In such a situation, 
the current is driven by a much more complicated mecha- 
nism (the coupling of reservoirs to the system) than sim- 
ply an external force, which can be nicely written as a 
term within the Hamiltonian of the system. Neverthe- 



less, the correlation function is ad hoc transferred to the 
density driven scenario by replacing the electrical current 
by the energy or heat current^. However, the justifica- 
tion of this replacement remains a conceptual problem 
here. 

One big advantage of this widely used approach is cer- 
tainly its computability after having diagonalized the sys- 
tem's Hamiltonian. A nice overview of results from the 
Kubo formula for spin models can be found, e.g., in the 
work by Heidrich-Meisner— (for further reading we sug- 
gest the comprehensive literature cited therein). How- 
ever, in most cases, a direct analytical solution for an in- 
finite system is not feasible and the interpretation of the 
results for finite systems seems to be not straight forward. 
For a finite system, the frequency dependent transport 
coefficient consists of numerous delta peaks with differ- 
ent weights at frequencies lo and is zero elsewhere. How 
to extract the dc-conductivity (interesting especially for 
the energy transport) of the finite system from this result 
or extrapolate the conductivity for the infinite one is a 
difficult questior>ii^i±&±2. 

A different approach to investigate the transport be- 
havior of a system is more connected to the experimen- 
tal measurement of heat conductivities: The system is 
directly coupled to heat baths of different temperatures 
within the theory of open quantum system a 18 i 19 ' 20 i 21 i 22 . 
That means that the Liouville von Neumann equation, 
describing the time evolution of the density operator of 
the system, is extended by incoherent damping terms 
simulating the influence of the heat baths. How to set up 
the correct dynamical equation here is highly nontrivial 
and involves the combination of many subtle approxi- 
mation schemes. In the case of an improper approach, 
the derivation can lead to mathematically correct, but 
physically irrelevant dynamical equations as discussed 
recently 2 ^. Having derived a proper quantum master 
equation (QME), the interpretation of the results for fi- 
nite systems is relatively easy: After finding the station- 
ary state of the dissipative dynamics, all interesting quan- 
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tities as currents and energy profiles are simply accessible 
by computing the expectation value of the respective op- 
erator. However, also this approach is restricted to finite 
systems since a complete analytical solution for larger 
systems is not available. Thus, the extrapolation to infi- 
nite systems needs a careful discussion to exclude errors 
due to the finite size of the investigated models as well. 

In the present paper, we will consider several model 
systems according to their transport properties. This is 
mainly done by the bath coupling method as discussed 
above, and by comparing with results from the Kubo for- 
mula. Let us start in Sect. [IT] with an introduction to the 
QME, the necessary observables and the Monte Carlo 
wave-function techniqu e) 24 ' 25 which is used to integrate 
the QME. Afterwards, we will present the results for sev- 
eral model systems: for chainlike systems in Sect. Mil and 
more complex ones in Sect. II Vl followed by our summary 
and conclusion. 



II. BACKGROUND 

A. Model system 

The considered system consists of N weakly coupled 
subunits described by the Hamiltonian 

N N-l 

H = H loc + H int = J2 +JY, hi ^ +1) ■ (!) 
/j— 1 /i— 1 

The first part H\ oc of the Hamiltonian contains the local 
spectra of the subunits. The second part H- m t describes 
the interaction between adjacent sites with the coupling 
strength J. Here, we require that the interaction is weak 
in the sense that the energy contained in the local part is 
much larger than the energy contained in the interaction 
part, (flioc) > (Him). 

More concretely, we will investigate one-dimensional 
(ID) or quasi-lD chains of two-level atoms or spin-1/2 
particles. Both, two-level atoms and spins, are described 
by the same algebra, and thus, it is convenient to use 
the Pauli operators {1, a x , a y , <r z } as a suitable operator 
basis here. The above mentioned weak coupling claim 
is fullfillcd by introducing a local Zeeman splitting ^ a z , 
where we require that f2 is much larger than the cou- 
pling constant J. Note that this weak internal coupling 
constraint is a necessary precondition for the validity of 
the master equation introduced below, i.e., we are not 
able to consider systems with Q approaching the same 
magnitude as J here. Hence, the models described by 
the Hamiltonian given in Eq. |T]) are not to be confused 
with spin chains in cuprates, e.g., where the local field 
is always close to zero, and thus, small compared to the 
coupling strength. However, the discussed models can be 
seen as spin chains in strong external fields, or simply as 
weakly coupled chains of two-level atoms as frequently 
considered in quantum optics and quantum information 
theory. 



To investigate the transport properties of these sys- 
tems, they will be explicitly coupled to independent en- 
vironments of different temperatures. Let us discuss the 
appropriate (QME)^i to describe this situation in the 
following Section. 

B. Lindblad Quantum Master Equation 

In general, the derivation of the QME from a mi- 
croscopic model 2,4 (a system coupled to an infinitely 
large environment) relies on some well known approxi- 
mation schemes the Marko v 26 ' 27 assumption, the Born 
approximation and the secular approximatio n 26 ' 28 . Re- 
cently, there was a discussion on how to derive a suitable 
Lindbla d 29 ' 30 QME in a nonequilibrium scenario^ 3 -, i.e., 
an equation to investigate transport in weakly coupled 
quantum systems. The Lindblad form of a QME defines 
a trace and hermiticity preserving, completely positive 
dynamical map 24 ' 31 , which thus retains all properties of 
the density operator at all times. In order to approach 
this dynamical equation the approximations are carefully 
carried out in a minimally invasive manner, to retain the 
central nonequilibrium properties of the model. It was 
shown that this nonequilibrium Lindblad QME is in very 
good accord with the results of the Redfield master- 
equation (non-Lindbladian) , contrary to the standard 
Lindblad QME in the weak coupling limit 2 ^. 

In a nonequilibrium investigation, one needs two heat 
baths at different temperatures locally coupled to the sys- 
tem, i.e., the heat baths couple only to a subunit at the 
edge of the system. The QME of such a situation yields 

^ = -i[H,p]+V L (p)+V R (p), (2) 

where the dissipator T>l refers to the left heat bath and 
Vr to the right one, depending on the full density oper- 
ator p of the system, i.e., the state of the chain described 
by the Hamiltonian ([1]). Both dissipators depend on the 
coupling strength A between system and bath as well as 
the temperature of the bath, respectively. Besides those 
incoherent damping terms, Eq. ^ contains a coherent 
part containing the Hamiltonian ([1]) of the system. 

The dissipator describing the heat bath coupled to a 
subunit at the edge of the system yields 

2 - 
V F (p) = £ ( 7 f)« (F k pF? - -[F}F k ,p]+) (3) 

k,l=l 

with F = L for the left and F — R for the right heat 
bath. The Lindblad operators Fk are given by 

U = ® i( 2 ) »...<» iW , (4a) 
L 2 = <7 W ® i( 2 > »...<» 1 w , (4b) 
4 = i (1) ®---®i (JV - 1) ®4 JV) , (4c) 

4 = i (1) ®---®i (iV - 1) ®^ JV) , (4d) 
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with the creation and annihilation operators a±. Here, 
the operators given in Eqs. (|4"aj) and (|4b|) belong to the 
left bath and those in Eqs. (|4c| and (|4djl to the right 
one. The coefficient matrices depend on the respective 
bath temperature (3p and are defined as 



If 



r F (fi) 



y/r F (n)v F {-si) 



N /r F (n)r F (-n) r F (-n) 



according to the rates 

r F (f>) 



All 



(5) 



(G) 



with the bath coupling strength A. This concrete form of 
the 7 F -matrices follows from a phenomenological ansatz 
for the spectral density of the environment, here chosen 
to be of Ohmi o 24 i 32 kind. 

A remarkable property of Eq. is that it can be 
brought into Lindblad form by diagonalizing the coef- 
ficient matrices jp. The complete dissipative part of 
Eq. ([2]) then reads 



V(p) 



4 

E 

fe=i 



pi- 



rn 



with Ek being linear combinations of the operators 
defined in Eqs. (|4all - (14djl and otk being non-negative num- 
bers. 

That it is indeed possible to derive a Lindbladian QME 
is very important here, since a standard stochastic unrav- 
elling of this special type of equation is feasible. Although 
extended stochastic schemes exist for the solution of gen- 
eral QMEs such as, e.g., the Redfield equatio n 33 i 34 i 35 , 
these methods have turned out to be less efficient, in 
general, than the standard approach. 



C. Observables and Fourier's Law 

The most interesting state of a nonequilibrium scenario 
is the local equilibrium state, i.e., the stationary state of 
the QME @. This state can be characterized by two 
central observables — the energy gradient and the energy 
current. Let us use 



fc(*0 = Tr{/jWp(£)} 



(8) 



as a local energy density at site p with pit) being the 
state of the system at time t. Since we are investigating 
internally weakly coupled subunits in the limit Q 3> J the 
local energy density is approximated by the local Hamil- 
tonian here. Therefore, we neglect completely the contri- 
butions to the local energy by the interaction. However, 
due to the smallness of J, these parts would be very small 
contributions to the above given energy density, and thus, 
would not dramatically change the results. 



In order to obtain a current operator between two ad- 
jacent sites in the system, we consider the time evolu- 
tion of the local energy operator given by the Heisenberg 
equation of motion for operators at site p 



4.hM=i[H,hM] + %-hM 
at at 



(9) 



Since is not explicitly time dependent the last term 
vanishes. Inserting Eq. ^ into ([9]) yields 



dt 



ftO*) = *([ftO'-i./O j £(/*)] + [h^ +1 \h^]\ . (10) 



Assuming that the local energy is a conserved quantity 
which is justified when Q J a discretized version of 
the continuity equation yields 

^-U^ = divJ = _ . (11) 

dt v ' 

By comparing Eq. (fl"0|) and Eq. (fTTj) we find for the cur- 
rent operator 



(12) 



Finally, the total energy current flowing from site p to 
site p + 1 is defined as 



(13) 



The celebrated Fourier's law (here a discrete version) 
states that in a proper diffusive situation, the current 
inside the system is proportional to the gradient, i.e., 



(14) 



here written in terms of energy current and energy gra- 
dient. If both current and gradient are equal at all 
sites p, and furthermore, the gradient is finite, a bulk 
conductivity^^ follows from 



(15) 



This is called normal or diffusive transport. On the other 
hand, if the gradient vanishes k diverges and the trans- 
port is called ballistic. However, that does not mean that 
the current diverges as well. Due to the resistivity at the 
contact to the heat bath (in our approach A) , the current 
will always remain finite. 

Even if we directly get a result in terms of normal or 
ballistic behavior for all finite systems here, a nonzero 
gradient in the finite system is not sufficient to deduce 
normal transport in the infinite one, too. The influence 
of the contact could dominate the investigation or long 
ballistic waves could be suppressed in the finite system. 
Thus, in order to obtain statements on the properties of 
the infinite system (bulk properties), it is important to 
investigate scaling properties as well. For normal trans- 
port behavior, both gradient and current must tend to 
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zero for infinitely large systems. Then and only then the 
system shows diffusive behavior. A finite current within 
an infinite system will always indicate ballistic transport 
behavior. 

Note that the current operator discussed in Eq. (fT2|) 
is essentially the standard spin current operator, mul- 
tiplied by the (large) Zeeman splitting f2. Due to the 
Zeeman splitting, there is an energy flow associated with 
any nonvanishing spin current. Only this energy flow is 
described by Eq. (fl"2"j) . i.e., it does not contain any en- 
ergy current that would be present even if Q was zero. 
Eventually, the stationary state of the QME @ will fea- 
ture such a nonvanishing spin current. Even though 
the z-component of the magnetization is conserved on 
the chain, the reservoirs may create and annihilate mag- 
netization in z-direction. Hence, we do not compen- 
sate for this current by applying an adequate magnetic 
gradient in the sense of Onsage r 36 i 37 (magnetic Seebeck 
effect 3 ^). Doing so, for large SI, the energy flow described 
by Eq. (fT2")l is the dominating part of the full energy cur- 
rent. 



D. Monte Carlo wave- function simulation 

In order to investigate the transport according to the 
QME requires the stationary solution p of Eq. From 
p all gradients and currents can be computed with Eq. ([5]) 
and Eq. (fT5)) . Unfortunately, Eq. (T2J) is an n 2 dimensional 
system of linear differential equations if n is the dimen- 
sion of the Hilbert space. To find the stationary state 
of this equation one has to diagonalize a n 2 x n 2 matrix 
which is restricted by the available memory. 

A very powerful technique to find the stationary state 
without diagonalizing the Liouvillian is based on the 
stochastic unraveling 24 of the QME. The basic idea is to 
depart from a statistical treatment by means of density 
operators and turn to a description in terms of stochas- 
tic wave functions. In fact, any QME of Lindblad form 
can equivalently be formulated in terms of a stochastic 
Schrodingcr equation (SSE) for the wave function \tp(t)) 



Am) 



iG(\m))m))dt 
' E k \m) 



iiw(*)>ii 



m)) dn k , (16) 



which describes a piecewise deterministic process in 
Hilbert space. The first term on the right hand side of 
Eq. (fTB|) describes the deterministic evolution generated 
by the nonlinear operator 

G{\m)) = h cS +^J2 <*k\\E k m))\\ 2 (17) 
A,- 

where we have introduced the non-Hermitian, effective 
Hamiltonian 



H eS = H - a k E\E k . 



(18) 



The second term in Eq. (JT6J) contains the Poisson incre- 
ments dn k £ {0, 1} which obey the following statistical 
properties 



(dnk) 
dn k dni 



\\E k \iP(t))\\ 2 dt 
S k i dn k . 



(19) 
(20) 



The stochastic process, defined by the SSE (fit))) , can 
be conveniently simulated by the following prescription. 
Starting from a normalized state, the first step of the 
unraveling procedure is to integrate the time-dependent 
Schrodinger equation according to the effective Hamil- 
tonian defined in Eq. (|18|) . Since it is not Hermi- 
tian the normalization of the state \^j(t)) decreases until 
{ip(t)\ip(t)} — rj 1 with r\ being a random number drawn 
from a uniform distribution on the interval {0, 1} at the 
beginning of the step. Subsequently, a jump k takes place 
according to the probability 



Pk = 



a fc ||W(*))|| 2 

E k *k\\E k \mw 



(21) 



Having identified the jump k, the state |^>(t)) is replaced 
by the normalized state 



\m) 



E k \m) 
\Eu\m)\\ 



(22) 



Afterwards, the algorithm starts from the beginning 
again with a deterministic evolution step. This proce- 
dure leads to one realization r of the stochastic process. 
Averaging over R — > oo realizations the time evolution of 
Eq. © is reproduced 2 ^. 

The expectation value of an observable A at time t can 
be estimated through 



Tr 



{ip(t)}«-^5> r (*)|i|^(«)> (23) 



in a finite ensemble of R realizations to arbitrary preci- 
sion. This is of huge practical importance, as one deals 
with wave functions with 0(n) elements instead of den- 
sity operators with 0(n 2 ) elements. Furthermore, if one 
is interested in the stationary state, ensemble averages 
can be replaced by time average a 39 ' 40 and one single re- 
alization suffices to determine the stationary expectation 
value 



1 T 

Tr{Ap} »A r = — J2(V(tk)\A\r(tk)) 



(24) 



fc=i 



with t k = to + kAt. It turns out that introducing this 
uniform time discretization, and allowing for jumps to 
occur at multiples of At only has several technical ad- 
vantages. However, one has to bear in mind, that this 
introduces an error of order O(At)^. A further prob- 
lem is that for At — > 0, the total number of timesteps T 
have to be increased in order to retain a sufficient num- 
ber of jumps in the average (|24p . This overall increase 
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FIG. 1: Chain of two- level atoms or spin- 1/2 particles coupled 
to heat baths of different temperature. 



of accuracy will be purchased at the cost of tedious com- 
putations. Nevertheless, in practice, the time-averaging 
procedure proves highly efficient, and results of sufficient 
accuracy could always be produced. It is further advis- 
able to discard the initial time evolution in the average in 
order to obtain reliable results, i.e., by choosing to ^> 0. 

In order to gain the standard deviation as a measure 
for the statistical error as well one should compute the 
stationary expectation value of R realizations. Thus, its 
mean is obtained by 



R 



and the standard deviation of the average yields 



1 R 



Af 



(25) 



(26) 



These errors are influenced by the chosen sampling in- 
terval At the neglected steps at the beginning to and the 
total amount of time steps T being averaged over. For all 
numerical results below we have chosen the parameters 
At = 1, < = 10 4 and T between 10 5 and 10 6 . For those 
settings the errors are surprisingly small already. 



III. CHAIN OF TWO LEVEL ATOMS 

First, we consider a chain of two-level atoms or spin- 
1/2 particles as depicted in Fig.[TJ In this case the local 
part of the Hamiltonian is just given by the mentioned 
Zeeman splitting of the individual spin 



9k 

2 



(27) 



with a splitting f2 M which may differ from site to site. 
Apart from that has to be large compared to the cou- 
pling constant J to remain in the weak coupling limit. 
The subunits are coupled by a generalized Heisenberg 
interaction 



0.2 0.4 0.6 0.8 1 

system number/TV 

FIG. 2: Local energies in a Heisenberg chain with length N = 
12 — 16. The system number is normalized by the chain length. 
The fit is carried out for chain length N — 16 excluding site 
one and 16. System parameters: J — 0.01, A = 1, Q = 1, 
A = 0.01, L = 0.5, p R = 0.25. 



For A^l the chain is called anisotropic chain and for 
A = the present model is equivalent to the XY-model 
(Forster coupling). In this case, by plugging Eq. (f27|) 
and the interaction given by Eq. (|28p into Eq. (fl2"|). the 
current operator yields 



M 5.(^+1) 



(29) 



(28) 



The above given system is coupled to heat baths of 
different temperatures. The left bath is set to the inverse 
temperature (3l = 0.5, and the hotter one at the right 
hand side is at fin — 0.25. Both baths couple with the 
same coupling strength A = 0.01 to the system. 

Having computed the stationary state of Eq. ([2]) by us- 
ing the method presented in Sect. HlDl one can compute 
both the stationary energy profile within the system and 
the current flowing through the system. In Fig. [2] the 
internal gradient is shown for an isotropic chain A = 1 
of TV = 12 — 16 spins according to the same constant 
local field Q = 1 and coupling strength J = 0.01. To 
show that the gradient is equivalent for the different sys- 
tem sizes we have normalized the chain length in Fig. [2] 
to one. The fit (line in Fig. \2§ is carried out for system 
size N = 16 excluding the sites one and 16, because of 
strong influences of the contacts. Even if the fit is done 
for system size 16 exclusively, all chains show the same 
gradient. However, the energy difference between adja- 
cent sites decreases for growing system sizes. The change 
in the internal gradient is shown in the upper diagram of 
Fig. [3] Here the gradient is plotted over the reciprocal 
chain length. The error bars refer to an average over the 
energy differences in all adjacent pairs of sites, where the 
first and last pairs have been neglected again as already 
done in Fig. O As can be seen from Fig. [3] the gradient 
decreases for larger systems until it approaches zero for 
an infinite chain which is in accordance with the expected 
behavior in the thermodynamical limit. 

The lower part of Fig. [3] shows the scaling behavior of 
the current through the system. Here, error bars refer 
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FIG. 3: Scaling properties of the Heisenberg chain N = 5— 16. 
Lines are fits carried out for chain length N = 6 — 16. System 
parameters: J = 0.01, A = 1, ft = 1, A = 0.01, /3l = 0.5, 
Pr = 0.25. 
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FIG. 4: Local energies in a XY-chain with length TV = 9 — 12. 
The system number is normalized by the chain length. The fit 
is carried out for chain length N = 12 excluding cite one and 
12. System parameters: J = 0.01, A = 0, O = 1, A = 0.01, 
/3 L = 0.5, (3 R = 0.25. 



to the failure produced by the stochastic algorithm given 
by the square root of Eq. (|26|) . The current decreases 
similar to the gradient, however, the extrapolation for 
the infinitely long chain does not approach zero. A finite 
current for an infinite system is a typical characteristic 
for ballistic transport behavior. According to the data 
shown in Fig. [3] one could eventually conclude finding 
ballistic transport in the Heisenberg chain. 

Figure [4] shows the local energy profile within the XY- 
model (A = 0). In comparison to Fig.[2]the profile within 
the systems is flat. According to Fourier's law (fT4|) this 
could be interpreted as ballistic behavior in the investi- 
gated finite models of different lengths (cf. discussion in 
Sect. Ill Cp . The local current between site /j, and fj, + 1 
remains finite although the gradient within the system 
vanishes. This local current is constant for all investi- 
gated system sizes and we find for the chosen parameters 
(5.33 ± 0.05) • 10 -4 . The results concerning the Heisen- 
berg chain and the XY-model are in accordance with 
some earlier results (for smaller systems) based on the 



FIG. 5: Dependence of the extrapolated value for the current 
of an infinitely long Heisenberg chain on the bath coupling 
strength A. 



full diagonalization of the Liouvillia n 18 ' . 

In Fig. [5] we investigate the dependence of the extrap- 
olated value of the current and the energy gradient of an 
infinitely long chain on the coupling strength A at the 
contact. In order to get comparable data and errors all 
other parameters are kept constant. A decrease of the ex- 
ternal coupling strength is combined with a decrease in 
the global decay time of the system and a drastic change 
of the jump probabilities (|21[) as well. To gain a proper 
expectation value from Eq. (|24[) with a rather small error 
it is crucial that the sampling time-step At of the contin- 
uous stochastic trajectory is chosen in a way that a suit- 
able amount of both coherent dynamics and stochastic 
jumps enter the average. That means, if At is too large, 
so that after each coherent step follows a jump already, 
the result of Eq. (|24|) will deteriorate. Thus, changing the 
external coupling strength would also require an adaption 
of the sampling parameter At. Furthermore, in case of 
a larger decay time of the system also the parameter to 
(initial neglect of data points) has to be increased. Thus, 
having fixed these parameters to get comparable data we 
are restricted to a small change of the external coupling 
strength only. For the finite system an increase in the ex- 
ternal coupling strength A denotes that a larger current 
is injected into the system, as follows from Fig. The 
resistance of the contact is decreased. Finally, this also 
results into a larger gradient within the system. Nev- 
ertheless, Fig. [5] shows that even if the results for finite 
systems changes drastically (especially for very small sys- 
tem sizes) the extrapolation for the infinite chain remains 
the same within the accuracy of the fit. 

Figure [5] shows the scaling behavior of the current for 
different values of the anisotropy A and Fig. [7] shows the 
scaling of the gradients, respectively. From the linear 
fits in Fig. [6] one could extrapolate the current within an 
infinitely long chain. This current is shown in Fig. [5] 
with dependence on A (A = 1 refers to the Heisenberg 
chain and A = to the XY-model) . Near the anisotropy 
A = 1.6 the current within the infinite system seems to 
vanish (cf. Fig. 0, i.e., the analysis at hand indicates 
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FIG. 6: Scaling behavior of the current in anisotropic Heisen- 
berg chains. System parameters: J = 0.01, O = 1, A = 0.01, 
Pl = 0.5, Pr = 0.25. 
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FIG. 7: Scaling behavior of the gradient in anisotropic Heisen- 
berg chains. System parameters: J = 0.01, O = 1, A = 0.01, 
Pl = 0.5, Pr = 0.25. 



normal transport behavior. Whether this is obtained for 
increasing A as well cannot be decided clearly from this 
analysis, but it seems to be probable that it remains dif- 
fusive for higher values of A. 

Having found diffusive behavior according to the Kubo 
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FIG. 8: Extrapolated current for the infinitely long chain with 
dependence on the anisotropy A. 
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FIG. 9: Local energies in an alternating local field Heisenberg 
chain N = 12 — 14. The system number is normalized by the 
chain length. The fit is carried out for chain length N = 14 
excluding the first and the last site. System parameters: J — 
0.01, A = 1, e = 0.02, A = 0.01, p L = 0.5, p R = 0.25. 



formula, it seems, nevertheless, unclear how to ex- 
tract the dc-conductivity from the behavior of the fi- 
nite system. Contrary to the Kubo investigation, a dc- 
conductivity directly follows from Eq. (|15p , in the present 
analysis, if we assume for the moment that the linear scal- 
ing of current and gradient found in Fig. [S]and[7|is also 
valid for larger systems. According to the small errors 
found in the above investigation, this assumption seems 
to be plausible. Thus, we are able to compute the con- 
ductivity of the infinite system for A = 1.6 using Eq. (fT5|) 
by dividing the slope of the current by the slope of the 
energy gradient directly finding = [2.34 ±0.08] • 10~ 2 . 
Here, the error follows from the uncertainty of the linear 
regression which is weighted already by the errors of the 
data points. 

Unfortunately, the models which can be investigated 
according to the suggested method are also restricted in 
size. The main restriction here is not the size of mem- 
ory, but the time one accepts to wait for the data. For 
the present technique, the computing time scales expo- 
nentially with the system size. Thus, investigations how 
disorder (random offset in the local field, random cou- 
plings) would change the above results are not available 
yet. 

Let us discuss another Heisenberg coupled chain of two 
level systems in the following. However, we will analyze 
an alternating local field in the following given by fl^ — 
1 + (— l) M e with fi = 2, 3, . . . , N— 1, i.e., no change in the 
field at the edge of the system. This keeps the contact 
unchanged even if e is varied. The gradient of such a 
system is depicted in Fig. [9] and shows a modulation in 
the local energies according to the change in the field. 
Nevertheless we have fitted the gradient by a least square 
fit, shown for N = 14 in Fig. [5J and used such fits to 
obtain data for the gradient in systems with different size 
(upper part of Fig. fTO]) . Even if the local profile is not as 
flat as before the current between adjacent sites is always 
the same. The scaling behavior of current and gradient 
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FIG. 10: Scaling properties of the alternating local field 
Heisenberg chain N = 5 — 14. Lines are fits carried out for 
chain length iV = 6 — 14. System parameters: J = 0.01, 
A = 1, e = 0.02, A = 0.01, j3 L = 0.5, (3 R = 0.25. 
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FIG. 11: Current of the infinite alternating local field Heisen- 
berg chain versus the parameter e characterising the change 
in the local field from site to site. 



TABLE I: Conductivity for the alternating local field chain 
in the diffusive regime. 
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1.29 ±0.04 
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FIG. 12: Ladder of two-level atoms or spins. Natural subunits 
are pairs of spins with an internal coupling strength J'. The 
coupling between the subunits is now given as two bonds in 
horizontal direction. 



ing two Heisenberg coupled spins with coupling strength 
J'. Thus, the local part of the Hamiltonian for subunit 
ll is described by 



1 + 1 



j'a 



(30) 



with the spin vector <r = {<j x , a y , a z }. Again we consider 
two-level systems with an energy splitting f2, here. The 
interaction between two adjacent sites ll is given by 



^n+i) _ ^ h3i<g>i(g><7i®i + i<g><7i®i( 



(31) 



is depicted in Fig. [TU] for e ~ 0.02. Again the current 
within the infinite chain seems to vanish, which could be 
a hint for normal transport. A scan over the parameter e 
is depicted in Fig. [TTJ As can be seen the current for the 
infinitely long chain roughly approaches zero at e = 0.02. 
For larger e it remains approximately zero within the 
stated accuracy. Thus, this investigation points towards 
normal transport behavior above e = 0.02. 

Using again Eq. (jT5|) and dividing the slope of the cur- 
rent by the slope of the gradient depicted, e.g., in Fig. [TO] 
we find the conductivities of the infinite system above 
e = 0.02 given in Tab. Q] Of course this is again only 
correct, if the scaling behavior remains the same as al- 
ready found in this finite size analysis for larger systems 
as well. 



IV. LADDER OF TWO LEVEL ATOMS 

The second class of models is a spin ladder introduced 
in Fig. 1121 In order to consider the transport in the 
model, the system is partitioned into subunits /i contain- 



Because of the weak internal couplings one may use a 
special type of bath contact called private bath here. This 
means that each spin at the edges of the system is coupled 
to its own "private" heat bath with temperature Pl at 
the left hand side and (3r at the right hand side. This 
is different from a more general approach where the two 
systems at the edge are viewed as one four level system 
each, to which the respective heat bath couples. The 
concept of private baths is only valid in the weak coupling 
limit again, i.e., J, J' <C fi. 

In Fig. [TO] we show the internal gradient of a ladder 
with N = 5 — 7 rungs, according to the same coupling 
strength in horizontal and vertical direction J — J' = 
0.01. Again we have normalized the length of the chain 
to fit all different chain lengths into the figure. As can 
be seen from the figure the system features a nice linear 
energy gradient, again. The scaling behavior of current 
and gradient for J = J' = 0.01 is depicted in Fig. [TO] 
There is no possibility to compute a reasonable error for 
the system sizes N ~ 3 and N = 4. This is due to the fact 
that tor N — 3 there is no pair without a contact to a heat 
bath and for TV = 4, there is just one energy difference 
in the center of the system away from the bath coupling. 
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FIG. 13: Local energies in a Heisenberg coupled ladder with 
N = 5 — 7 rungs. The system number is normalized by the 
number of rungs. The fit is carried out for N — 7 excluding 
site 1 and 7. System parameters: J = 0.01, J' = 0.01, Q = 1, 
A = 0.01, p L = 0.5, /3r = 0.25. 



FIG. 15: Current of the infinite Heisenberg ladder versus the 
vertical coupling strength J'. 



but physically not interpretable. 



o 

•3 



0.0008 - data 
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reciprocal chain length 

FIG. 14: Scaling properties of the Heisenberg ladder with 
N = 2 — 7 rungs. Lines are fits carried out for rungs N = 5 — 7. 
System parameters: J = 0.01, J' = 0.01, Q = 1, A = 0.01, 
P L = 0.5, p R = 0.25. 



Nevertheless, we have plotted those system sizes within 
Fig. 1141 The errors for larger system sizes are surprisingly 
small, because of the large total system size. Due to the 
too small available system size and the strong influence of 
the baths at the edges, the gradient in the small systems 
is essentially different from the gradient in larger ones. 
However, the current in small systems is already close to 
the linear fit for larger system sizes. 

Figure 1 1 51 shows the current for the infinite system ex- 
tracted of the scaling analysis in dependence of the verti- 
cal coupling strength J'. According to this investigation 
we do not find normal transport in any of the considered 
ladder models. The coupling strength J' = J already 
seems to be the minimum. For J' = we get two com- 
pletely independent Heisenberg chains with a maximum 
length of seven spins. Larger values than J' = 0.02 could 
be evaluated, however, results would be questionable be- 
cause the weak coupling limit is violated. Since the weak 
coupling limit is crucial for the derivation of the underly- 
ing QME, the results would be mathematically correct, 



V. CONCLUSION 

In the present paper we have studied several models 
of two- level atoms or spin- 1/2 particles coupled to heat 
baths of different temperatures. The investigation has 
been based on a recently derived quantum master equa- 
tion (QME) which simulates the noncquilibrium situation 
properly. This QME is of Lindblad form and can, thus, 
be efficiently unravelled using a standard Monte Carlo 
wave-function technique. The significant advantage of 
such an approach is its applicability to larger systems 
in comparison to the restricted system sizes which can 
be investigated by a direct diagonalization of the Liou- 
villian. This follows from the fact that the stochastic 
unraveling deals with wave functions rather than density 
operators. All interesting quantities are, here, given as 
mean values over stochastic trajectories and can be eval- 
uated to arbitrary precision by adjusting the amount of 
timesteps averaged over. Here, we are interested mainly 
in the energy profile within the system and the energy 
current through it. The analysis of both current and 
gradient in dependence of the system size gives informa- 
tion on the underlying transport behavior, which could 
be, in principle, of ballistic or normal diffusive nature. 

For finite systems, the method at hand always leads to 
an easily interpretable result, in terms of currents, energy 
profiles, and the resulting conductivities. This is mainly 
a result of the concrete design of the method by a direct 
contact of the system with heat baths. The extrapola- 
tion to infinite system sizes, to extract bulk properties 
from the analysis of finite systems, is done by a careful 
scaling analysis of both gradients and currents. Finally, 
this provides an indication of the type of transport in an 
infinite probe of the model system as well. 

We have analyzed a multitude of different concrete spin 
models here. Those systems consist of weakly coupled 
two-level atoms or spins. In the first part, we have in- 
vestigated chains according to different coupling models 
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TABLE II: Comparison between results for the energy trans- 
port (as defined by Eq. p2p ) within different model systems 
featuring a finite local splitting f2 (finite magnetic field), as 
obtained from the Kubo formula and the bath coupling ap- 
proach. Results for the Kubo formula are taken from Refs. 13 
and 43. 



Infinite Model 


Bath 


Koo[10- 2 ] 


Kubo 


generalized Heisenberg 

XY-model, A = 


ball. 




ball. 


Heisenberg, A = 1 


ball. 




ball. 


A < 1.6 


ball. 




ball. 


A = 1.6 


diff. 


(2.34 ±0.08) 


ball. 


A > 1.6 


diff. 




ball. 


Alternating chain 

e < 0.02 


ball. 




??? 


e = 0.02 


diff. 


(1.29 ±0.04) 


??? 


e > 0.02 


diff. 




??? 


Heisenberg ladder 


ball. 


??? 



and local fields. The consideration was centered around 
the generalized Heisenberg chain, i.e., a Heisenberg chain 
with different ZZ-coupling strengths according to the pa- 
rameter A within the model. Among those, one promi- 
nent model is the XY-model with vanishing A featuring 
ballistic transport in both the finite as well as the infinite 
model. Despite the relatively small system sizes investi- 
gated here, there are some hints for normal transport in 
the models as well: The scaling analysis of the anisotropic 
model with A = 1.6 features a vanishing current for in- 
finite chain size. Such a finding is typically connected 
to diffusive behavior. Second, the Heisenberg chain with 
alternating local field shows a zero current above some 
threshold dependent on the difference of the local fields. 

Besides those investigations concerning the transport 
behavior, we have also considered the dependence of the 
results on the bath contact. Here, we find that the ex- 



trapolated results do not crucially depend on the respec- 
tive coupling strength. 

In the second part we have analysed a Heisenberg lad- 
der with different coupling strength in horizontal and ver- 
tical directions. According to this investigation of very 
short ladders, no evidence for normal diffusive behavior 
with dependence on the coupling strength in vertical di- 
rection is found. The extrapolated current remains finite 
over the complete accessible parameter space. 

In a nutshell, the comparison between our results and 
results based on the Kubo theory can eventually be sum- 
marized in Tab. [TTJ 

The heat bath coupling approach to transport behavior 
presented in this paper features some significant advan- 
tages: Besides the direct determination of current and 
gradient in the finite system, and thus, of the conductiv- 
ity in any concrete finite situation, the method also al- 
lows us to extract the conductivity of the infinite model 
by means of an extrapolation. This extrapolation relies 
on the linear scaling behavior of current and gradient to- 
ward zero, which is extracted from the analysis of finite 
system sizes. Thus, the direct coupling of reservoirs to 
model systems substantially improves the understanding 
of the transport behavior of such models on both small 
and infinite scales. 
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